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General relativistic simulations for the merger of binary neutron stars are performed as an exten- 
sion of a previous work [l[. We prepare binary neutron stars with a large initial orbital separation 
and employ the moving-puncture formulation, which enables to follow merger and ringdown phases 
for a long time, even after black hole formation. For modeling inspiraling neutron stars, which 
should be composed of cold neutron stars, the Akmal-Pandharipande-Ravenhall (APR) equation of 
state (EOS) is adopted. After the onset of merger, the hybrid-type EOS is used; i.e., the cold and 
thermal parts are given by the APR and F-law EOSs, respectively. Three equal-mass binaries, each 
with mass I.AMq, 1.45M0, and 1.5M0, and two unequal-mass binaries with mass, 1.3 vs l.GM© and 
1.35 vs 1.65M0, are prepared. We focus primarily on the black hole formation case, and explore 
mass and spin of the black hole, mass of disks which surround the black hole, and gravitational 
waves emitted during the black hole formation. We find that (i) the black hole is promptly formed 
if total mass of the system initially satisfies mo > 2.QMq\ (ii) for the systems of mo = 2.9-3.OM0 
and of mass ratio ~ 0.8, the mass of disks which surround the formed black hole is 0. 006-0. 02M0 ; 
(iii) the spin of the formed black hole is 0.78 ± 0.02 when a black hole is formed after the merger 
in the dynamical time scale. This value depends weakly on the total mass and mass ratio, and is 
about 0.1 larger than that of a black hole formed from nonspinning binary black holes; (iv) for the 
black-hole formation case, Fourier spectrum shape of gravitational waves emitted in the merger and 
ringdown phases has a universal qualitative feature irrespective of the total mass and mass ratio, 
but quantitatively, the spectrum reflects the parameters of the binary neutron stars. 

PACS numbers: 04.25.D-, 04.30.-w, 04.40.Dg 



I. INTRODUCTION 



Coalescence of binary neutron stars is one of the most promising sources for kilometer-size laser interferometric 
detectors such as the LIGO [3, i], GEO VIRGO 0, Q, and TAMA 7]. Latest statistical estimate indicates 
that detection rate of gravitational waves from binary neutron stars will be 1 event per ~ 40-300 years for the first- 
generation interferometric detectors and ~ 10-100 events per year for the advanced detectors This suggests 
that gravitational waves from binary neutron stars will be detected within the next decade. 

The merger of the binary neutron stars has also been proposed as a likely candidate for the central engine of short 
7-ray bursts (GRBs) [l^jllli- The observational facts that short GRBs have a cosmological origin (e.g., Ref. [T^ ) 
indicate that the central engine supplies a large amount of energy > 10^* ergs in a very short time scale ^ 0.1-1 
s [isj . According to a standard scenario based on the merger hypothesis, a stellar-mass black hole surrounded by a 
hot and massive disk (or torus) should be formed after the merger. Possible relevant processes to extract the energy 
of this black hole-accretion disk system for launching a relativistic jet are neutrino-anti neutrino annihilation and 
magnetically driven mechanisms. 

Recent semi-analytic calculations (e.g., Ref. [3]) show that an accretion rate of M > O.IM0/S is required for neu- 
trino annihilation models to achieve sufficiently high energy efficiency. Also, recent numerical studies (e.g., Ref. (Tsj ) 
suggest that if the disk had a mass > O.OIMq, it could supply the required energy by neutrino radiation. Recent 
general relativistic magnetohydrodynamic simulations jlJI] of (Kerr) black hole-accretion disk system find that the 
rotational energy of the black hole can be extracted via magnetohydrodynamics processes as a form of the Poynting 
fiux. In particular the results are in general agreement with the predictions of Blandford and Znajek [iTj . In this case, 
the black hole is required to be rapidly rotating for efficient energy extraction. This fact motivates to calculate final 
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mass and spin of the black hole formed after the merger, and to clarify whether such a massive disk can be formed, 
and to clarify what type of the progenitors are necessary for formation of a system composed of a black hole and a 
massive disk. 

For theoretically studying the late inspiral, merger, and ringdown phases of the binary neutron stars, numerical 
relativity is the unique approach. Until quite recently, there has been no longterm general relativistic simulation 
that self-consistently clarifies the inspiral and merger phases because of limitation of the computational resources 
or difficult y in simulating a black-hole spacetime, although a number of simulations have been done for qualitative 
studies 

0, in, ^20, 22,, 23,, 24, 25]. The most crucial drawbacks in the previous works are summarized as follows 
(but see Refs. [32, ISJI); (i) the simulations were not able to be continued for a long time after formation of a black 
hole and/or (ii) the simulations were short-term for the inspiral phase; the inspiral motion of the binary neutron 
stars is followed only for ~ 1-2 orbits. The simulations were usually started with a quasiequilibrium state which is 
obtained assuming that approaching velocity between two neutron stars is zero, that is not realistic. As a result of this 
treatment, non-zero approaching velocity at the onset of merger is not correctly taken into account, and moreover, 
effects of non-zero eccentricity could play an unfavored role |34l |. 

In the present work, we perform an improved simulation overcoming these drawbacks; (i) we adopt a moving- 
puncture approach [2a . [27j . which enables us to evolve black hole spacetimes for an arbitrarily long time; (ii) we 
prepare binary neutron stars in quasiequilibrium states of a large separation as the initial condition. In the chosen 
initial data, the binary neutron stars spend ~ 4 orbits in the inspiral phase before the onset of merger, and hence, 
approximately correct non-zero approaching velocity and nearly zero eccentricity results. Furthermore, we add a non- 
zero approaching velocity at < = 0. The magnitude of the approaching velocity can be estimated by a post-Newtonian 
(PN) analysis (see Sec. |IT] in detail). Adding the ap pro aching velocity suppresses an artificial orbital eccentricity 
caused by an incompleteness of the initial conditions [28|. In addition, a non- uniform grid with a sufficiently large 
computational domain [26l . [29l [sol . [3l| is employed to perform longterm accurate simulations with a relatively low 
computational cost. 

A longterm simulation of binary neutron stars, in which the inspiral phase is followed for 3-5 orbits, has very 
recently been performed by Baiotti et al. using the polytropic and F-law EOSs for modeling the neutron stars [33j . 




With these simple EOSs, they self-consistently investigated the inspiral, merger, and ringdown processes. However, as 
is well known, such EOSs are oversimplified for modeling the neutron stars, and hence, the results are only qualitative. 
In this paper, we perform a longterm simulation not with such simple EOSs but by adopting a nuclear-theory-based 
EOS. We focus in particular on quantitatively clarifying the formation process of a black hole for the case that it is 
formed promptly (i.e., in the dynamical time scale ~ 1-2 ms) after the onset of merger. More specifically, the primary 
purpose of this paper is (1) to determine the condition for prompt black hole formation, (2) to determine the final 
mass and spin of the black hole formed after the merger, (3) to clarify quantitative features of gravitational waves 
emitted in the merger and ringdown phases, and (4) to estimate the mass of disks surrounding the formed black hole. 

Here, some words of explanation are necessary about nuclear-theory-based EOS (often referred to as realistic EOS). 
There are two types of nuclear-theory-based EOSs for modeling the neutron stars. One is zero-temperature nuclear 
EOS such as the Akmal-Pandharipande-Ravenhall (APR) EOS ^55]. This type of EOSs could be suitable for modeling 
not so young neutron stars, such as neutron stars in binary neutron stars just before the merger, for which the thermal 
energy per nucleon is much smaller than the Fermi energy. The other is finite-temperature EOS such as the Shen's 
EOS This type of EOSs should be employed for studying high-energy phenomena such as supernova core collapse 
and the merger phase of binary neutron stars. 

One problem faced when choosing EOSs is that no one knows the truly realistic EOS for high-density nuclear matter. 
This implies that for studying the inspiral and merger phases of binary neutron stars by numerical simulation, it is 
necessary to systematically perform many simulations adopting a variety of nuclear-theory-based EOSs. Indeed, one 
of the important roles in gravitational-wave detection is to constrain nuclear EOS from detected gravitational waves 
0. For this purpose, numerical simulation for a variety of EOSs has to be performed for deriving all the possible 
gravitational waveforms. There are many cold EOSs that have been proposed [37| . and a systematic theoretical 
survey of the gravitational waveforms for clarifying their dependence on the EOSs is possible. By contrast, there are 
only a few finite-temperature EOSs (sgI. [ssj. This implies that a systematic study is not possible when adopting the 
finite-temperature EOSs. 

In the series of our papers P,[2l],[22|, we have found that shock heating during the merger phase increases the thermal 
energy of the neutron stars. This indicates that adopting finite-temperature EOSs is desirable for a sophisticated 
study of the merger process. However, the shock heating is not very strong during the merger and the thermal energy 
per nucleon after the shock heating is only ~ 10% of the Fermi energy of the nuclear matter in the central region of 
the merged object (see Sec. IV). Thus, the thermal energy is not essential at least for short-term evolution of the 
main body of the merged object. In particular, the thermal energy plays a minor role in determining the criterion for 
prompt black hole formation at the merger and its process. Because of this reason, in the previous two papers [U. [23j. 
we have adopted a hybrid EOS in which the pressure and thermal energy are divided into cold and hot parts. The 
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cold part is determined by the zero-temperature nuclear-theory-based EOSs, and the hot part is simply modeled by 
a F-law EOS. In this method, the main part (cold part) of the EOS is appropriately modeled, although a minor part 
(finite-temperature part) is modeled in an approximate manner. The merit of this EOS is that a variety of the cold 
EOSs can be systematically employed for the simulation. The hybrid EOS may not be appropriate for studying the 
longterm evolution of a hypermassive neutron star which is formed when the total mass of the system is not large 
enough for prompt black hole formation. The reason for this is that in the hypermassive neutron star, longterm shock 
heating and cooling via neutrino emission are likely to play a role. However, for studying prompt black hole formation, 
which is the main subject of this paper, the hybrid EOS seems to be acceptable because the roles of the shock heating 
and neutrino cooling do not seem to be essential. Thus, following previous papers [H, [1^, numerical simulations are 
performed, employing the hybrid EOS. Specifically, we adopt the APR EOS for the cold part following Ref. 

This paper is organized as follows. In Sec.|lTl we describe the formalism for numerical solution of the Einstein and of 
relativistic-hydrodynamics equations, the EOS adopted in this paper, and the method for extraction of gravitational 
waves. The initial condition for binary neutron stars and grid structure for the simulation are summarized in Sec. IIIII 
In Sec. IIV| numerical results are presented, focusing in particular on the case that a black hole is formed in the 
dynamical time scale after the onset of merger. Sec. |V] is devoted to a summary. Throughout this paper, unless 
otherwise stated, we adopt the geometrical units in which G = c = 1 where G and c are the gravitational constant 
and the speed of light. Greek and Latin indices denote the spacetime and space components, respectively. 



II. FORMULATION 
A. Numerical methods 

Our formulation for the fully general relativistic simulation is the same as in Refs. [H, [H, [13], to which the reader 
may refer for details of the basic equations. 

For solving the Einstein evolution equations, we use the original version of the Baumgarte-Shapiro-Shibata- 
Nakamura (BSSN) formalism [s^: We evolve the conformal factor, cj) = (ln7)/12, the trace part of the extrinsic 
curvature, K, the conformal three-metric, jij = 7^^^'^7y , the tracefree extrinsic curvature. Ay = ^^^^^{Kij—Kjij 
and a three-auxiliary variable, Fj = S^^dj^ik- Here 7ij is the three-metric, Kij the extrinsic curvature, 7 = det(7y ), 
and K = Kijj^^ . As in Ref. [30|, we evolve the conformal factor 0, not the inverse of ■0, because the cell-centered 
grid is adopted in our code, and hence, the coordinate singularity at the puncture is avoided in moving puncture 
frameworks [l^, [l^] • 

For the conditions of the lapse, a, and the shift vector, we adopt a dynamical gauge condition in the following 
forms, 

{dt-P'd,)\iia = -2K, (2.1) 
dtl3' = 0.75f^Fj + MdtFj), (2.2) 

where Ai denotes the time step in the numerical simulations, and the second term on the right-hand side of Eq. (j2.2p 
is introduced for stabilizing the numerical computations. The gauge condition (|2.2p . which was originally proposed 
in Ref. [7 0|, i s slightly different from that usually used in the moving puncture framework (see, e.g., [1^). We note 
that Ref. [33] shows that this gauge is as suitable as popular one for simulating black hole spacetimes. 

The numerical scheme for solving the Einstein equations is essentially the same as that in Ref. [1^. We use the 
fourth-order finite difference scheme in the spatial direction and a fourth-order Runge-Kutta scheme in the time 
integration, where the advection terms such as (3^di4> are evaluated by a fourth-order upwind scheme, as proposed in 
Ref. [27] . In a previous uni-grid simulation, we use a third-order scheme for the time integration [s^ . We have found 
that the fourth-order scheme can give more accurate results, and hence, updated the scheme for this work. 

The location and properties of the black hole, such as the area and the circumferential radii, are determined by 
analyzing an apparent horizon. Our method for finding the apparent horizon is described in Refs. [40ll4]|. From the 
area, and polar and equatorial circumferential proper lengths, we infer the mass and spin of the black hole. 

The numerical code for the hydrodynamics is the same as that in Refs. [2^, |33]: As the variables to be evolved, we 
adopt = pau^e^'^ , Ui = hui, and e* = hau^ — Pj (pau*), where p is the rest-mass density, Ui is the three-component 
of the four velocity, m* is the time component of the four velocity, P is the pressure, h is the specific enthalpy defined 
by /i = 1 -|- e -|- P/p, and e is the specific internal energy. To handle advection terms in the hydrodynamic equations, 
a high-resolution central scheme [43] is adopted with a third-order piecewise parabolic interpolation and with a steep 
min-mod limiter. In the present work, the limiter parameter, 6, is set to be 2.5 (see Ref. [i^] for detail about our 
interpolation scheme and the parameter 5). 
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B. Equations of state 

Following Refs. P, H^l) we adopt a hybrid EOS for modeling the EOS of neutron stars. In this EOS, the pressure 
and the specific internal energy are written in the form 

P = Pcold + Pth, (2.3) 

£ = ecoid + eth, (2.4) 

where Pcoid and EcoM are the cold (zero-temperature) parts, and are written as functions of rest-mass density p. In 
general, any nuclear-theory-based EOS for zero-temperature nuclear matter can be employed for assigning PcoM and 
Scold- In this paper, we adopt the APR EOS [s^, for which P and e are tabulated as functions of the baryon rest-mass 
density for a wide density range. We use a fitting formula for the data in the range, 10"'^'' g/cm'^ < P ^ 10^^ g/cm"^. 
The method of the fitting was first developed in Ref. [i^l and slightly modified in Ref. [J, [2^ . We adopt the fitting 
parameters listed in Table I of Ref. [l| . 

Pth and £th in Eqs. ()2.3|) and \2A\ denote thermal (finite-temperature) parts which are zero in the absence of 
shocks, but become non-zero if shocks are formed. Specifically, they have finite values after the merger sets in. During 
the simulation, p and e are determined from the evolved variables p*, e*, and the normalization relation of the four- 
velocity. Then, £th is determined by e — £coid (p) , and subsequently, the thermal part of the pressure Pth is related to 
the specific thermal energy eth = e — Ecoid as 

Pth = (Eth - l)p£th, (2.5) 

where Fth is an adiabatic constant for which we set Pth = 2 taking into account the fact that the EOSs for high-density 
nuclear matter are stiff. 



C. Extracting Gravitational Waves 



To extract gravitational waves from numerical data, we compute the outgoing component of the Newman-Penrose 
quantity \E'4 (e.g., Refs. [13, IH, for detail). From ^'4, loss rates of energy and angular momentum carried by 
gravitational waves are computed by 



dE_ 



= lim 



r 

16^ 



(i(cos 9)dip 



t 2' 



—— — lim Re 

dt r^QO 



r 



—— i> d(cos9)dip 



d^-^idt' 




^idt'dt" 



(2.6) 
(2.7) 



where § d(cos 9)dip denotes an integral on two surface of a constant coordinate radius and #4 is the complex conjugate 
of ^'4. In numerical simulation, the surface integral is performed for several radii near the outer boundaries, and we 
check that the resulting gravitational waveforms depend only weakly on the extracted radii. The outer boundaries 
along each axis are located dX r — L ^ \q where Aq is wave length of gravitational waves emitted by inspiraling binary 
neutron stars at i = and gravitational waves are extracted for r « 0.7-0.95L (cf. Table [iTl maximum extraction 
radii are denotes by r^^ in this table). 

The amplitude and phase of gravitational waves extracted at several radii should be extrapolated to obtain those 
at infinity for deriving precise waveforms, as often done in the simulation for binary black holes [28l. ItH. WX FtsI . 
However, numerical waveforms computed in hydrodynamic simulation are not as precise as those in the simulations 
for vacuum spacetime, because the order of the accuracy is reduced at shocks, discontinuities, and places where the 
gradient of hydrodynamic quantities is large in the standard numerical hydrodynamics. This implies that numerical 
error is primarily determined by such sources, and hence, the finiteness of the extraction radius does not become 
the main source of the numerical error. To illustrate this fact, we generate Fig. [1] which shows the time evolution of 
gravitational wave amplitude for a typical model in this paper. The amplitude A(t) is defined from ^4 as 



r*4 



A{t)e 



(2.8) 



where r is the extraction radius and </> the phase. This plot shows that the amplitude depends very weakly on the 
extraction radius; relative difference among three results is much smaller than 1%, typical error size induced by the 
poor resolution of hydrodynamics mentioned above. Thus, we do not use extrapolation of amplitude and wave phase 
in this work. 
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For calculating the two polarization modes /i+,x from ^"4, we perform the time integration of twice with 
appropriate choice of integration constants and subtraction of unphysical drift, perhaps associated primarily with the 
drift of the mass center of the system, that often occurs in hydrodynamic simulation due to accumulation of numerical 
error. Specifically, whenever the time integration is performed, we subtract a function of the form 02^^ + ait + qq 
where 09-02 denote constants which are determined by the least-square fitting to the original numerical data. 

From a time sequence of dE/dt and dJ/dt, we compute total radiated energy and angular momentum by time 
integration as 

AE = Jdt^, (2.9) 
AJ=[dt!^. (2.10) 



D. Comparison with PN approximation 



Gravitational waveforms in the inspiral phase computed in numerical simulation should be compared with those by 
the PN approximation, because it provides an accurate waveform if the orbital velocity v is not extremely relativistic 
(i.e., V < 0.3c ). Assuming that binary components are point masses and their quasicircular orbits evolve adiabatically, 
orbital evolution of the binary stars can be analytically determined with 3PN accuracy and gravitational waveforms 
are derived with the 3.5PN accuracy [Tsj. Recently, several phenomenological prescription for improving the PN 
results, which is applicable even to highly relativistic orbits of u ~ 0.3c, have been proposed. Among them, high- 
accuracy simulations for equal-mass and nonspinning binary black holes [2^| have proven that the so-called Taylor T4 
formula provides the orbital evolution and gravitational waveforms with a high accuracy at least up to about one orbit 
before the onset of the merger. The Taylor T4 formula appears to be a good approximation also for unequal-mass 
and nonspinning binaries [46| (but see Ref. [811] for a depression of the accuracy for the Taylor T4 in spinning cases). 
Because we focus only on nonspinning binaries, it is safe to assume that the Taylor T4 formula would be a good 
approximate formula, and thus, we explore a match between the numerical and Taylor-T4's waveforms for the inspiral 
phase. 

One issue in comparison in the present context is that the effect of tidal deformation of neutron stars, which is not 
taken into account in the Taylor T4 formula, plays an important role for close orbits. This implies that the numerical 
waveforms should not agree with the waveforms by the Taylor T4 formula for such orbits. Thus, we compare two 
waveforms for < 0.04 (see Sec. Ill A for definition of mo). 

More specifically, in comparison, two waveforms have to appropriately align as done in Refs. [H, |7l|, |7§|. Our 
procedure in this work is as follows: First we calculate the orbital angular velocity by = d^/dt. Then, we 
determine the reference time at which two waveforms give the same value of il. Recent studies in the context of the 
binary black holes have shown that the reference value of fi should be < 0.1/mo because beyond this value the PN 
approximation breaks down [H, IH, |7l|, ■ For the binary neutron stars, the value of 0. 1/ mo is even too large because 
the compactness of neutron stars is much smaller than that of black holes and hence the merger already started at 
such a high value. Also, the effect of tidal deformation of neutron stars plays an important role for mo^l > 0.04. Due 
to this reason, the reference value of is chosen to be moil « 0.04 for all the runs in this paper. 



III. INITIAL MODEL AND SIMULATION SETTING 



A. initial model 



Except for the orbit just before the merger, binary neutron stars are in a quasicircular orbit because the time scale 
of gravitational radiation reaction at Newtonian order ~ (5/64)17"^ (Mf?)^/^ (see, e.g., (60j ') is several times longer 
than the orbital period. Hence, following our previous works P, [12, lH, [13, [H, [4l[, we adopt binary neutron stars 
in quasiequilibrium states as initial conditions. The quasiequilibrium states are computed in the so-called conformal- 
flatness formalism for the Einstein equations [i^l . The irrotational velocity field is assumed because it is considered 
to be a realistic velocity field for coalescing binar y n eutron stars in nature [3 [1^ • We employ numerical solutions 
computed by a code in the LORENE library [sol Tsil. [s^. [ssj . Table |T] lists the several key quantities for the models 
adopted in this paper. We select three equal-mass (APR1414, APR145145, and APR1515) and two unequal-mass 
models (APR1316 and APR135165). Note that the numbers in the model name denote the ADM mass of the two 
neutron stars in isolation (e.g., for APR1316, mass of two neutron stars in isolation is mi = 1.3AfQ and m2 = 1.6AIq). 



In the following, we use mo, Mq, and as the sum of masses of two neutron stars in isolation (mp = mi + m2), 
initial total ADM mass, and total rest mass of the system, respectively. 



B. grid setting 

In the simulation, the cell-centered Cartesian, (x,y,z), grid is adopted. In these coordinates, we can avoid the 
situation that the location of the puncture (which always stays on the z — plane) coincides with one of the grid 
points. Equatorial plane symmetry is also assumed. The computational domain of — L < x < L, —L < y < i, and 
< z < L is covered by the grid size (2iV, 27V, N) for {x, y, z), where L and iV are constants. Following Refs. [29l[30j. 
we adopt a nonuniform grid as follows; an inner domain is covered with a uniform grid of spacing Ax and with 
the grid size, {2No,2No, Nq). Outside this inner domain, the grid spacing is increased according to the relation, 
^tanh[(i — A'o)/Az]Ax, where i denotes the i-th grid point in each positive direction, and Nq, Ai, and f are constants. 
Then, the location of i-th grid, x''{i), in each direction is 

fc,., r (* + l/2)Aa: < i < No , . 

^ '~\{i + l/2)Ax + ^AiAxlog[cosh{{i--NQ)/Ai}] i > Nq ^"^'^^ 

and x^{—i — 1) = —x^{i), where i = 0, 1, • • • iV for x^ = x, y, and z. The chosen parameters of the grid structure for 
each simulation are listed in Table |TT1 

For investigating convergence of numerical results, we perform simulations with three grid resolutions for all the 
models; labels "L" , "M" , and "H" denote the low, medium, and high grid resolutions. We note that for all the runs, 
the grid resolution around the neutron stars is much better than that in the previous work in which the major 
diameter of the neutron stars is covered only by 45 grid points. Since it is covered by « 60-80 grid points, we expect 
that numerical results in this paper are much more accurate than those in the previous work. In physical units, the 
grid spacing around the neutron stars are about 200 meter for the highest grid resolution. On the other hand, the 
grid spacing in a wave zone, where gravitational waves are extracted, is ~ 4 km. Radius of apparent horizon of black 
hole finally formed is covered by about ten grid points. 

For high-resolution run (run "H"), about 200 GBytes computational memory is necessary and 240 CPU hours was 
spent using 512 processors on Cray XT4 system at the Center for Computational Astrophysics (CfCA) in the National 
Astronomical Observatory of Japan (NAOJ). For run "L", the computational time is ~ 100 CPU hours using the 
same processors. 



C. approaching velocity 

In computing quasiequilibrium states in the conformal-flatness formalism, approaching velocity, which should be 
present in reality due to gravitational radiation reaction, is not taken into account. Lack of the approaching velocity 
induces a non-zero orbital eccentricity and resulting modulation in gravitational waveforms (e.g., Ref. [2^). If a 
simulation is started from an initial condition in which the initial orbital separation is sufRciently large, the value 
of the approaching velocity is negligible, and hence, its lack is not a serious problem; in such case, the approaching 
velocity settles to a correct one and eccentricity approaches zero in a few orbits, because of gravitational radiation 
reaction [H^. However, in the present choice of the initial condition, the initial orbital separation is not sufRciently 
large. Thus, we initially add an approaching velocity to improve the quality of the initial condition. For calculating an 
approximate value of the approaching velocity, we assume that the two neutron stars may be approximated by point 
particles and the so-called Taylor T4 formula 28|^45j is used for predicting their orbital motion. In the following, we 



describe our method (see relevant works [45. l75l l76j in simulations of binary black holes): Assuming that the density 
maxima of the neutron stars are located along x-axis at i = 0, the approaching velocity may be calculated by 

X _ dri2 _ GniQ ^712 , , 

~ dt ~ 7?2c2 dt ' ^-^-^^ 

where 712 = GrriQ / ri2C^ and ri2 is a coordinate orbital separation in the harmonic gauge, mg is the total mass defined 
by the sum of masses of two neutron stars at a state when they are in isolation. Here, we recover G and c to explicitly 
clarify that 712 is dimensionless. 712 is written by a gauge-invariant quantity, x = (Gmof2/c'^)^/'^, at 3PN order as 
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TABLE I: List of several key quantities for the initial data of binary neutron stars in quasiequilibrium 
state. The ADM mass of each star when they are in isolation (mi and 7712), the maximum baryon 
rest-mass density for each star, the baryon mass ratio Qm = M«2/M«i, the total baryon rest mass 
M«, the total ADM mass Mq, nondimensional spin parameter Jq/Mq, orbital period Pq, and the 
orbital angular velocity in units of M^^ , MoQ,o, where f2o denotes initial orbital angular velocity. 



Model 


mi, m2(MQ) 


p (10^5 g/cm-') 


Qm 


M.(M0) 




Jo /MS 


Po(ms) 


Mono 


APR1414 


1.40, 1.40 


0.887, 0.887 


1.00 


3.106 


2.771 


0.983 


3.185 


0.0269 


APR145145 


1.45, 1.45 


0.935, 0.935 


1.00 


3.232 


2.870 


0.983 


3.299 


0.0269 


APR1515 


1.50, 1.50 


0.961, 0.961 


1.00 


3.359 


2.969 


0.983 


3.412 


0.0269 


APR1316 


1.30, 1.60 


0.864, 1.015 


0.7943 


3.238 


2.870 


0.976 


3.299 


0.0269 


APR135165 


1.35, 1.65 


0.887, 1.045 


0.7992 


3.365 


2.970 


0.978 


3.412 


0.0269 



where J7 denotes the orbital frequency, v = mim2/mg is the ratio of the reduced mass to the total mass with mi (7712) 
being the mass of the star 1 (2) in isolation, and rg is "logarithmic barycenter" [18]. In the Taylor T4 formula [28| , 
the evolution equation of x is derived in an adiabatic approximation as 



dx 



16c3 
5Gmo 



487 
168 



X + 47ra;' 



3/2 



274229 
72576 ^ 



254 



^5/2 



178384023737 1475 



3353011200 



192 



1712 
T05 



■lE 



856, ; 

^ln(16.) 



3310 
189"' 



,7/2 



(3.4) 



where 7e is Euler's constant. From Eqs. p.2p - p.4p . the approaching velocity is derived. 

We note that this approaching velocity is not gauge-invariant and has a strict meaning only in the harmonic gauge 
condition. Because the initial condition is not obtained in this gauge condition, it would be necessary to carry out 
a coordinate transformation to obtain the approaching velocity in our chosen gauge condition. However, the initial 
condition is obtained in the conformal flatness formalism, and thus, the gauge condition is not clear, because the 
metric components are oversimplified in this formalism. Thus, in this work, we use the approaching velocity defined 
in Eq. (j3.2p with no modification. After calculating w^, the four velocity is modified from its quasiequilibrium value 

u^^'^'' by adding a correction u'^"" — we^'^v^ /a with w = au* = (1 -I- "f^^u^'^'^^u'^''^^)^^'^ where we drop a correction in a 
and due to the change of Ui, because they only give higher-order corrections. Then, we uniformly add the specific 
momentum hu'^°" to each neutron star and reimpose the Hamiltonian and momentum constraints. Note that it is 
nontrivial to relate velocity to momentum in general relativity (but see Ref. f77l| for discussion about this issue). 

Figure [2] plots the evolution of a coordinate orbital separation /^/^ with/without approaching velocity for model 
APR1515. Here, / is defined by 



1 = 



(3.5) 



where ly- = / pau^x'^x^ yFf(Px with u* being a time component of the four velocity. In the absence of the approaching 
velocity, / does not decrease monotonically but oscillates with time (the dotted curve in Fig. [2|). In the presence of 
the approaching velocity, this oscillation is suppressed (the solid curve in Fig. This figure illustrates the advantage 
of our treatment. It should be noted that the orbital eccentricity is not completely suppressed even in this method. 
As discussed in Sec. IIV C 1[ indeed, gravitational waveforms in an early phase slightly disagree with that calculated 
in the Taylor T4 framework (e.g., angular velocity computed by gravitational waves does not increase monotonically 
but has a modulation; see Sec. IIV C ip . However, the radiation reaction circularizes the binary orbit and gravitational 
waves in the late inspiral phase agree approximately with the prediction by the Taylor T4 formula in a better manner 
(cf. Fig. 13 (b)). 



IV. RESULT 



A. General feature for merger process 

We have already performed simulations for binary neutron stars employing nuclear-thcory-based EOSs [J, [23l |. 
Although the qualitative feature for the merger process found in the present work is the same as in the previous 
works, we here summarize generic feature of the merger again. 
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TABLE II: Parameters for the grid structure employed in the numerical simulation. The grid number 
for covering one positive direction (A'^), that for the inner uniform grid zone (Aq), the parameters 
for nonuniform-grid domain (Ai,^), the approximate grid number for covering the major diameter 
of massive neutron star (Lns)i the ratio of the outer grid spacing to the wavelength of fundamental 
quasinormal mode of the formed black hole (Aqnm), the ratio of the location of outer boundaries 
along each axis to the initial gravitational wave length (Ao = 7r/r2o) and maximum radius of 
gravitational wave extraction, rex, in units of kilometer (and in units of L). The last column shows 
the final outcome; BH and NS denote a black hole and neutron star, respectively. 



Model 


A 


Ao 


Ai 




Lj^s/Ax 


L/\o 






Outcome 


APR1414L 


224 


114 


30 


15 


60 


0.97 


— 


4.15E+2 


(0.90) 


NS 


APR1414M 


234 


120 


30 


17.5 


70 


0.98 


— 


4.35E+2 


(0.94) 


NS 


APR1414H 


253 


139 


30 


20 


80 


0.97 


— 


4.28E+2 


(0.93) 


NS 


APR145145L 


224 


110 


30 


20 


60 


1.24 


7.0 


4.41E+2 


(0.72) 


BH 


APR145145M 


252 


123 


30 


20.5 


70 


1.24 


8.3 


5.67E+2 


(0.92) 


BH 


APR145145H 


282 


141 


30 


21 


80 


1.24 


9.2 


5.80E+2 


(0.95) 


BH 


APR1515L 


224 


114 


30 


19.5 


60 


1.10 


7.3 


5.06E+2 


(0.90) 


BH 


APR1515M 


250 


128 


30 


20 


70 


1.09 


8.3 


5.24E+2 


(0.94) 


BH 


APR1515H 


280 


145 


30 


21 


80 


1.10 


9.6 


5.42E+2 


(0.94) 


BH 


APR1316L 


239 


130 


30 


19 


60 


1.11 


7.7 


5.05E+2 


(0.92) 


BH 


APR1316M 


278 


160 


30 


20 


70 


1.10 


8.6 


4.93E+2 


(0.91) 


BH 


APR1316H 


300 


170 


30 


21 


80 


1.11 


9.5 


5.13E+2 


(0.94) 


BH 


APR135165L 


240 


130 


30 


20 


60 


1.11 


7.8 


5.15E+2 


(0.90) 


BH 


APR135165M 


278 


155 


30 


20.5 


70 


1.11 


8.9 


5.14E+2 


(0.90) 


BH 


APR135165H 


310 


175 


30 


21 


80 


1.10 


9.9 


5.19E+2 


(0.92) 


BH 



Figure [3] plots the evolution of the minimum value of the lapse function, Ofmin, and maximum baryon rest-mass 
density, p^ax, for all the models studied in this paper. For models APR145145, APR1515, APR1316, and APR135165 
for which a black hole is formed in the dynamical time scale p^^/^), ckmin (Pmax) decreases (increases) monotonically 
after the onset of merger. In the high-resolution runs for these models, two neutron stars come into the first contact at 
t ~ 9 ms. This merger time is slightly underestimated because of the effect of finite grid resolution, but the numerical 
results in the chosen grid resolution is in a convergent regime as discussed in Sec. lIVCTI For all the cases, an apparent 
horizon is formed when amin reaches ~ 0.03. For model APR1414, ckmin (pmax) steeply decreases (increases) after the 
onset of merger, but then, they start oscillating and eventually settle down to relaxed values. This indicates that the 
outcomes is a hypermassive neutron star, for which the baryon rest mass is ~ 30% larger than the maximum allowed 
mass of the spherical neutron stars [ssj . 

Figures |3Hni display the evolution of density contour curves for the rest mass in the equatorial plane of the inspiral, 
merger, and ringdown phases for runs APR1414H, APR1515H, and APR1316H, respectively. For the case of the 
equal-mass binaries, two neutron stars are tidally deformed in a noticeable manner only just before the onset of 
merger (see Figs. |l]and[5])- If the total mass is not large enough to form a black hole in the dynamical time scale 
after the onset of the merger, a hypermassive neutron star of nonaxisymmetric structure is formed in the central 
region. In this case, spiral arms are formed in the outer region (see Fig. U]), and subsequently, they wind around the 
formed hypermassive neutron star, generating shocks in its outer region. Because of angular momentum dissipation 
by gravitational radiation, hydrodynamic interaction, and hydrodynamic transport process of angular momentum, 
the density and specific angular momentum of the hypermassive neutron star are redistributed during the subsequent 
evolution and eventually it relaxes to a moderately ellipsoidal configuration. Here, the ellipticity remains because it is 
rapidly rotating and also the EOS is stiff enough as mentioned in Refs. [lillll- a result, the hypermassive neutron 
star emits gravitational waves subsequently, and it secularly evolves due to gravitational radiation reaction. 

For the case when the total mass is large enough, the merged object collapses to a black hole in the dynamical time 
scale after the onset of merger (see Fig. [SJ. In the equal- mass case, nearly all the material are swallowed by the black 
hole in ~ 1 ms. The resulting final outcome is a black hole surrounded by a tiny disk of mass < lO^'^Af0. 

In the unequal-mass case, the less massive neutron star is tidally deformed ~ 1 orbit before the onset of merger (see 
the second panel of Fig. [6]). Then, mass shedding occurs, and as a result, the material of the less massive neutron star 
accrete onto the massive companion. During the merger, it is highly tidally deformed, and thus, an efhcient angular 
momentum transport occurs. Due to this, the material in the outer region of the less massive neutron star spread 
outward to form a spiral arm. This process helps formation of accretion disks around the formed black hole. The disk 
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will survive for a time scale much longer than the dynamical time scale as shown in Sec. IIVBI 

Figure [7] displays the density contour curves for the rest-mass density and the local velocity field on the x = 
plane for runs APR1414H, APR145145H, APR1316H, and APR135165H. These also show that (i) a hypermassive 
neutron star of ellipsoidal shape is the outcome for model APR1414 (panel (a)), (ii) a black hole with tiny surrounding 
material is the outcome for model APR145145 (panel (b)), and (iii) a black hole surrounded by disks is the outcome 
for models APR1316 and APR135165 (panels (c) and (d)). 

All these features are qualitatively the same as those reported in Ref. [1]. However, in the previous work, we 
employed initial conditions for which the initial separation is much smaller than that in the present work, as mor^o 
0.05 where JIq denotes initial orbital angular velocity. In the present work, mgrio ~ 0.027, and this modification 
changes the results quantitatively. Furthermore, we performed the simulations for different total masses, and this 
leads to a new quantitative finding. In the remaining part of this subsection, we summarize the updated results. 

As reported in the previous works [U, HI], the final outcome formed after the merger (a black hole or neutron star) 
is primarily determined by a relation between the initial total mass mg and the threshold mass Mthr and by the mass 
ratio of the binary for a given EOS: The binary neutron stars of mo > A/thr collapse to a black hole after the onset 
of merger in dynamical time scale ~ 1 ms. On the other hand, a hypermassive neutron star is formed for mo < Mthr, 
at least for a time much longer than the dynamical time scale. The threshold mass depends weakly on the mass ratio 
of the binary neutron stars. 

The threshold mass, Mthr, depends strongly on the EOSs (on the stiffness of the EOS), and stiffer EOSs give larger 
values of Afthr- We reported in the previous work that for the APR EOS, 2.8Mq < Mthr < ^-OMq. In the present 
work, we find that the value of Afthr is in a narrower range between 2.8Mq and 2.9M0 via the study for models 
APR1414, APR145145, and APR1316. (Note that for models APR1414, APR145145, and APR1515, mo = 2.8Mq, 
2.9M0, and S.OM©, respectively) 

In the previous paper we reported that the disk mass around the black hole is ~ 4 x 1O~^M0 for the equal-mass 
binary neutron star of total mass mo = 3M0 and ^ O.OO3M0 for the binary neutron star of mass 1.35 and 1.65M0. 
We find that the mass of the disk at i - Uh = 3 ms is w 7 x IO^^Mq for model APR1515H and « 6 x IO^^Mq for 
model APR135165 (see Fig. [TUl) . Here, ^ah is the time when the apparent horizon is first formed. Thus, the disk 
mass is corrected by a factor of several, although the order of magnitude does not change. Also, we did not clarify the 
dependence of the disk mass on the total mass mo in the previous work. For model APR1316, we find that the disk 
mass at i — ^ah = 3 ms is 2.4 x lO~^Af0, and thus, the disk mass depends strongly on the total mass of the system. 
More details about the merger process, final outcome, and disk mass will be discussed in subsequent subsections for 
each model separately. 



1. APRI414 

The total mass of this model is slightly smaller than the threshold mass, i.e., mo < Mthr, and hence, a compact 
hypermassive neutron star is formed. Because the total mass is close to Mthr in this model, the merged object first 
becomes very compact; e.g., amin and Pmax reach ~ 0.2 and ~ 2 x 10^^ g/cm^ soon after the first contact (see Fig. 
13] (a)). Hence, the self-gravity is only slightly insufficient for inducing collapse to a black hole. The merged object 
subsequently bounces back to become an oscillating hypermassive neutron star. The quasiradial oscillation is repeated 
for several times (see Fig. [3]), and then the oscillation is damped due to shock dissipation and the hypermassive neutron 
star relaxes to a quasisteady ellipsoidal configuration [1], as found from the last panel of Fig. [Hand Fig. [7] (a). These 
figures show that axial ratio of the two major axes on the equatorial plane is ~ 0.9, and that of the polar coordinate 
radius rp to the equatorial one To is T-p/vc ~ 0.7. Thus, the ellipticity is not negligible. Such high cUipticity seems to 
reflect the fact that the hypermassive neutron star is rapidly rotating [56j . 

Figure [5] (a) plots angular velocity profile of the hypermassive neutron star in the late phase along x and y axes. 
This shows that the hypermassive neutron star rotates rapidly and differentially. In the central region, the rotational 
period is ~ 0.5 ms, comparable to the dynamical time scale. This implies that centrifugal force around the central 
region plays an important role for supporting self-gravity. Note that the mass of the hypermassive neutron star is 
by « 30% larger than the maximum allowed mass of spherical neutron stars of the APR EOS (which is ~ 2.2M0), 
and hence in the absence of rotation, the hypermassive neutron star would collapse to a black hole. Thus, the rapid 
rotation near the central region seems to be an essential agent for supporting its strong self-gravity. (We note that 
thermal energy generated by shock heating during the merger phase in part plays a role for supporting the self-gravity; 
see discussion below.) 

As we will discuss in Sec. IIV CI the hypermassive neutron star continuously emits gravitational waves and loses 
angular momentum because it has a nonaxisymmetric shape and rapid rotation. This indicates that it may eventually 
collapse to a black hole after a substantial fraction of the angular momentum is dissipated from the central region. The 
lifetime of the hypermassive neutron star may be estimated by calculating the time scale for the angular momentum 
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loss due to the gravitational- wave emission. At the end of the simulation {t ~ 20 ms), the angular momentum and 
its dissipation rate are J ~ 0.74Jo ~ 4.9 x lO^^g cm^ s~^ and dJ/dt ~ 6.7 x 10"'^ g cm^s~^, respectively. The 
lifetime of the hypermassive neutron star may be estimated by J/{dJ/dt) ^ 700 ms. dJ/dt (namely the amplitude 
of gravitational waves) gradually decreases with time, and hence, this time scale should be regarded as the shortest 
one. However, the decrease time scale of dJ/dt is not as short as the dynamical time scale, and thus, the estimated 
time scale is likely to be correct within the factor of 2-3. If this estimation is correct, the hypermassive neutron 
star would collapse to a black hole in a few seconds (but see discussion below for other possibilities). Other physical 
processes, which are not taken into account in this work, could also contribute to dissipating and/or transporting 
angular momentum: Because the hypermassive neutron star rotates differentially as shown in Fig. [8l magnetic fields 
might be amplified by the magnetorotational instability and/or magnetic winding [stI. [ssI. [s^ . As a result, angular 
momentum may be transported efficiently, leading the hypermassive neutron star to collapse to a black hole. If the 
time scale for the magnetic processes is shorter than the emission time scale of gravitational waves, they would be 
the main agent for inducing gravitational collapse to a black hole. 

As discussed above, the hypermassive neutron star eventually collapses to a black hole in any scenario. Because 
the hypermassive neutron star has a spread envelope (see the last three panels of Fig. 3), the final black hole formed 
after the gravitational collapse may be surrounded by an accretion disk. To qualitatively estimate the outcome, we 
generate Fig. [8] (b) which shows the evolution of the mass spectrum as a function of the specific angular momentum 
M:f{j) where j = hu^p is the specific angular momentum of each fiuid element. Here, M^(j) is defined as an integrated 
baryon mass of fluid elements with j > j'; 



Figure [5] (b) shows that the value of j for most of the fluid elements is smaller than 2Mo for t <20 ms. However, the 
fraction of the mass element of j > Mq increases with time. This indicates that angular momentum is transported 
outward due to nonaxisymmetric hydrodynamic interaction. 

The ADM mass and angular momentum at the end of the simulation {t w 20 ms) are w 0.97Afo and 0.74Jo, 
respectively, and thus, the nondimensional spin parameter of the system is 0.77. Then, assume that the final mass 
and spin of the black hole would be Mbh ~ 0.97Mo and a ^ 0.77. Specific angular momentum at innermost stable 
circular orbit (ISCO), jisco, of such a rotating black hole is given by « 2.45Mo (e.g., Ref. [lO])- This suggests that 
a fluid element of j > 2.45Mo would be able to form a disk surrounding the black hole. Figure [5] (b) shows that any 
mass element does not have specific angular momentum large enough to form the disk. However, the profile of the 
mass spectrum quickly changes with time by the hydrodynamic angular momentum transport, as mentioned above. 
The estimated lifetime of the hypermassive neutron star is a few seconds, and much longer than 20 ms. This suggests 
that a substantial fraction of the fluid elements may form an accretion disk. 

Finally, we touch on thermal effects on the evolution of the hypermassive neutron star. Because shocks are generated 
during the merger and in the subsequent dynamical phase, the hypermassive neutron star is heated up and as a result, 
the specific thermal energy, Sth, becomes nonzero. To clarify the role of the generated thermal energy, we plot profiles 
of Eth, P, and Pth/-Pcoid along x and y axes in Fig. [5](c)~(e). Figure[5](c) shows that the value of £th/c^ is 0.002-0.05 for 
a region of p > p„uc ~ 2 x 10^* g/cm"^, where pnuc denotes the nuclear density, and for the region of subnuclear density, 
it is 0.01-0.02. Assuming that the matter field is composed of neutron gas and thermal radiation, the temperature of 
the region with p > 10^^ g/cni"^ (above which the optical depth for neutrino transport would be larger than unity and 
the cooling due to the neutrino emission would not be efficient) is approximately calculated to give a high temperature 
as 7.2 X 10^°(eth/0.01c^) K. Nevertheless, the thermal pressure Pth in the central region of p > 10^^ g/cm^ is only 
~ 1-2% of Pcoid as shown in Fig. 7 (e), and hence, the thermal pressure plays a minor role for supporting the 
self-gravity in the central region. This indicates that centrifugal force plays a more important role for supporting the 
self-gravity of the hypermassive neutron star in its early evolution phase. By contrast, for a region of p < Pnuc, the 
thermal pressure is larger than Pcoid- Namely, the effect of the thermal energy plays an important role for determining 
the profile in the envelope of the hypermassive neutron star. 

Although the ratio Pth/^coid is small in the central region, it is not a negligible value and the thermal pressure seems 
to contribute in part to supporting the self-gravity. This suggests that even after dissipation and/or transportation of 
angular momentum via gravitational radiation or magnetic-field effects, the hypermassive neutron star of relatively 
small mass may not collapse to a black hole because of the presence of the thermal energy. If so, the collapse will set 
in after the thermal energy is dissipated via neutrino cooling. The cooling time scale by neutrino emission from the 
central region of the hypermassive neutron star is likely to be of order 10-100 s (e.g., chapters 11 and 18 of Ref. 60]), 
and thus, the lifetime of the hypermassive neutron star of relatively low mass may be rather long. 

In all the above scenarios, the hypermassive neutron star eventually collapses to a black hole after a substantial 
fraction of angular momentum is dissipated. This indicates that the resulting black hole will not be rotating as rapidly 
as the black holes promptly formed after the onset of merger (cf. Sec. IIVBI) . 




(4.1) 
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For the envelop of the hypermassive neutron star, £th ~ 0.01-0.02 and the temperature of this region is high, ^ 10^^ 
K. Because of its high temperature and relatively low density, such region is subject to a large amount of thermal 
neutrino emission. In particular, for the region of p < lO^^g/cm'^, neutrinos are not trapped by the matter but escape 
freely, so cooling will proceed rapidly. Thus, the thermal energy decreases on a time scale, shorter than the dynamical 
time scale, until the cooling time scale becomes as long as the dynamical one. This point should be explored in the 
future, incorporating finite-temperature EOSs and the neutrino emission, although such work is beyond scope of this 
paper. 

2. APR145145,APR1515 

In these equal-mass massive models, a black hole is formed after the onset of merger in dynamical time scale ~ 1 ms, 
and most of the material falls into the black hole (e.g., the last two panels of Fig. [5]). Figure [TUl (a) plots the evolution 
of Mr>rAHi which denotes the rest mass of baryon located outside the apparent horizon, for runs APR145145II and 
APR1515H. Here, Mr^rAH is defined by 

Mr>rAH = / P*d^x, (4.2) 

where rAH[= ''ah(^, f)] denotes the radius of the apparent horizon. 

This figure shows that more than 99% of the fluid elements are swallowed by the black hole within ^0.1 ms 
after the formation of the apparent horizon. The primary reason for this rapid infalling is that the specific angular 
momentum jini at the onset of merger is too small for all the fluid elements: The maximum value of the specific 
angular momentum at the onset of merger, jmax, is about 1.25Mo for both the models. As discussed in Sec. IIVBI 
the final value of the black-hole spin is estimated to be « 0.8 for both the models. This predicts that the specific 
angular momentum at ISCO around the formed black hole is « 2.4Mo, which is much larger than jmax- Moreover, the 
angular momentum transport works inefficiently because of equal-mass symmetry (see Fig. [5]) and also the specific 
angular momentum decreases due to the gravitational-wave emission during the merger. All these facts indicate that 
formation of massive disks surrounding the black hole is unlikely. 

The disk mass for model APR1515 is smaller than that for model APR145145. The likely reason for this is that the 
object formed just after the onset of merger for model APR1515 is more compact than that for model APR145145. 
As a result, (i) the dynamical time scale becomes shorter and the time duration for which the angular momentum 
transport works does as well; (ii) dissipation of angular momentum by gravitational waves is larger and disk formation 
becomes less likely. This dependence of the disk mass on the total mass is also found in the merger of unequal-mass 
binary neutron stars (see Sec. IIV A3p . 

Accretion time scale for the disk and increase time scale of the area of the apparent horizon arc much longer than 
the dynamical time scale (i.e., rotational time scale of the disk) for t — ^ah ^ 1-5 ms, as shown in Fig. [TO] (b). 
Therefore, the final state formed for these models is a rapidly rotating black hole surrounded by a quasisteady disk 
of tiny mass (see Fig. [Zl(b)). 

We compare the result for the disk mass in the present and previous simulations for model APR1515. In Ref. fl|, we 
reported that the final disk mass is '--^ 4 x 10~'*Mo, whereas it is ~ 7 x 1O~^M0 for the present work. This difference 
seems to originate simply from the fact that we evaluated the disk mass soon after its formation {t — ^ah = 0.5 ms) 
in Ref. [l|; i.e., the disk mass is evaluated before the disk relaxes to a quasisteady state. In the present work, we can 
evaluate it at t — ^ah = 3.0 ms, at which the disk is in a quasisteady state, because we employ the moving puncture 
method which enables us to perform a longtcrm simulation even after black hole formation. 

3. APR1316,APR135165 

As shown in Ref. [T] and in Fig. 5, the less massive neutron star is significantly tidally deformed in the final inspiral 
phase of unequal-mass binary neutron stars, and then, the merger occurs. At the onset of merger, the less massive 
star is highly elongated and hence its outer part subsequently forms a large spiral arm. The merger in the central part 
simultaneously proceeds and, for the case that too > Mthr, a black hole is formed around the center in a dynamical 
time scale ~ 1 ms. One interesting feature in the early stage of the merger (before gravitational collapse to a black 
hole) is that material in the inner part of the less massive star slips through the surface of the companion neutron 
star and forms a small spiral arm. As a result, two asymmetric spiral arms are formed around the central object 
which subsequently collapses to a black hole (see the fourth panel of Fig. [S]). Due to the nonaxisymmetric structure, 
the angular momentum is transported outwards in the spiral arms. Because the rotation velocity of the small spiral 
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arm around the central object is faster than that of the large spiral arm, they collide within one orbit. As a result, 
shocks are formed and the material in the spiral arms is heated up. In this mechanism, kinetic energy of the material 
is converted into thermal energy, and then a part of the fluid elements, which lose the kinetic energy, are swallowed 
by the black hole (see Fig. [TU)) . The fluid elements which escape from falling into the black hole eventually form a 
disk around the black hole. All these features were already found in Ref. [1], but in the previous work, the simulation 
was not able to continue for a time long enough to determine the final state because the simulation crashed before a 
quasisteady state was reached. As Fig. [TU] (b) indicates, the present simulations, by contrast, have been done until 
the final state (composed of a black hole and quasisteady accretion disk) is reached, and the conclusive statement 
becomes possible. 

Figure [7| (c) and (d) also show that the final outcomes for models APR1316 and APR135165 are a rotating black 
hole surrounded by a disk. The disk mass for runs APR1316H and APR135165H evaluated at t — ^ah = 3.0 ms is 
about 2.4 X 10"^ and 6.4 x IQ-^Mg, respectively (see Fig. [TUKa) and Table [m|. The hkely reason for this significant 
difference is that for the less massive binary neutron star, the system at the onset of merger is less compact and has 
a longer dynamical time scale for transporting angular momentum outward more efficiently, as already mentioned in 
Sec. IIVA2I 

The disk mass computed for three grid resolutions are described in Table Hill We find that the disk mass does not 
systematically converge and magnitude of the error seems to be by 50%. The likely reason for the slow convergence 
is that the grid resolution in the outer domain (where the fluid elements spread) is not systematically improved as the 
grid resolution in the inner domain is improved and disk mass depends sensitively on spurious numerical transport 
process of angular momentum (see, e.g., Ref. [zl] for discussion about slow convergence of disk mass in black hole- 
neutron star merger). However, all the numerical results are derived in a fairly narrow range irrespective of the grid 
resolution, and the order of magnitude of the disk mass does not depend on the grid resolution. 

In the previous paper |l|, the merger of unequal-mass binary neutron stars collapsing to a black hole was studied 
only for toq = 3M0. We reported that for the mass ratio Qm ^ 0.8, the rest mass of the disk surrounding the black 
hole was likely to be smaller than O.OIM©. The present results indicate that the disk mass depends sensitively on the 
total mass of the system, and even for the case that the mass ratio is not much smaller than unity, a system composed 
of a black hole and massive disk may be an outcome after the merger, if the total mass is close to Afthr- 

For model APR1316, the disk mass is larger than O.OIA/0 with a high maximum rest- mass density ~ 10^^ g/cm'^ 
in its inner region, as shown in Fig. [5) Also, the thermal energy is generated by shock heating, resulting in a high 
specific internal energy, e/c^ ~ 0.01 which implies a high matter temperature as 7.2 x 10^'^(eth/0.01c^) K (see Fig. [S]). 
All these properties are favorable for producing a large amount of high-energy neutrinos, from which electron-positron 
pairs and gamma-rays may be generated jisj. Thus, the final outcome is a candidate for the central engine of short 
GRBs. 



B. Black hole mass and spin 

As mentioned in the previous subsections, the final outcomes after the merger for models APR145145, APR1515, 
APR1316, and APR135165 are a rotating black hole. We here determine the black-hole mass, AfsH.f, and spin, af, 
using the same methods as those used in Refs. (30ll32|. 

There are at least two methods for (approximately) estimating the black- hole mass and three methods for estimating 
the black-hole spin. In the first method for estimating the black-hole mass and spin, we use the conservation laws. 
The energy conservation law is approximately written as 

MBH,f = Mo ~ Mr>r^^ " A^, (4.3) 

where Mq is the initial ADM mass of the system. In Eq. (|4.3p . we neglect the binding energy between the black hole 
and surrounding matter, but it is likely to be at most 10% of Air>rAH and thus a minor correction. 
Conservation law of angular momentum is approximately written as 

^BH,f = Jo - .A>rAH - A J, (4.4) 

where JeH.f is angular momentum of the black hole and Jq is total angular momentum at f = 0. Jr>rAH approximately 
denotes angular momentum of the material located outside apparent horizon which is defined by 

Jr>rAH = / p^hu^d^x, (4.5) 

where is the (/^-component of the four velocity. We note that Jr>rA-a is strictly equal to the angular momentum of 
material only for stationary axisymmetric systems. Thus, this should be adopted only for the case that the system 
approximately relaxes to a stationary axisymmetric spacetime. 
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From AfBH,f and JsH.f , a nondimensional spin parameter is defined by 

Hereafter, we refer to this spin as Ofi. 

In the second method, the mass and spin of the black holes are determined using geometrical properties of apparent 
horizon. When the system relaxes to a quasisteady state, the black-hole mass may be approximately estimated from 
equatorial circumferential length Ce of apparent horizon, because Ce/47r is equal to Men.f for Kerr black holes. This 
value also gives an approximate value of the black-hole mass even in the presence of surrounding torus [6lj . 

Next, we assume that area of apparent horizon, Aah, obeys the same relation as that of Kerr black holes; 



= " 2 ■ ^ ^ 

Then, from ^ah and Ce, a black-hole spin may be estimated. Hereafter, we refer to this spin as af2. Note that this 
method should be used only in the case that the system relaxes to a quasisteady state, because Eq. (|4.7p holds only 
for such spacetime. 

A black- hole spin is also estimated from polar and equatorial circumferential radii of apparent horizon, Cp and Ce- 
For Kerr black holes, the ratio Cp/Ce is given by a known function composed only of the black-hole spin as 



Cp 




(4.8) 



where f+ = 1 + V 1 — . The black-hole spin determined from Cp/Ce is referred to as aa in the following. 

All the results for the black-hole mass and spin are summarized in Table Hill We find that Mbu,! and Ce/47r agree 
within 0.5% error for all the models. This indicates that both quantities at least approximately denote the black-hole 
mass and that we obtain the black-hole mass within ~ 0.5% error. 

The spin parameters af2 and Ofs agree within 0.01 for all the models. However, the spin parameters Ofi does 
not agree well with af2 and ots. The typical size of the difference is 0.06-0.07, 0.05-0.07, and 0.04-0.06 for low-, 
medium-, and high-grid resolutions, respectively. Because the magnitude of the difference decreases systematically 
with the improvement of the grid resolution, this discrepancy originates primarily from numerical error associated 
with the finite grid resolution. As discussed in Ref. ^32j . poor grid resolution enhances spurious dissipation of angular 
momentum and shortens the time duration of the inspiral phase. This leads to the underestimation of A J. Indeed, 
the time duration of the inspiral phase and AJ increase with improvement of the grid resolution, and as a result, 
the value of Ofi systematically decreases (see Table HI). By contrast, Cp/Ce and Aah depend weakly on the grid 
resolution, and so do af2 and aa. This suggests that the convergent value of the black- hole spin is close to af2 and 
at3, and therefore, we conclude that for all the models, the black-hole spin is af_~ 0.78 ± 0.02. 

In the merger of equal-mass, nonspinning binary black holes, at is « 0.69 [H, |4^, which is by about 0.1 smaller 
than that for the merger of binary neutron stars. This difference arises primarily from the magnitude of dJ/dt in 
the final phase of coalescence. In the merger of binary black holes, a significant fraction of angular momentum 
is dissipated by gravitational radiation during the last inspiral, merger, and ringdown phases, because a highly 
nonaxisymmetric state is accompanied by a highly compact state from the last orbit due to the high compactness 
of the black holes. By contrast, compactness of the neutron stars is a factor of ~ 5-7 smaller, and as a result, 
such a highly nonaxisymmetric and compact state is not achieved for the binary neutron stars. Indeed, the angular 
momentum loss rate by gravitational waves during the last phases is much smaller than that of the binary black holes. 

Finally, we comment on dependence of the black hole mass and spin on the approaching velocity. As discussed in 
Sec. IllICi our prescription of the approaching velocity can reduce the initial eccentricity. By comparing the result for 
run APR1515L with/ without approaching velocity, we find that the spin and mass of the final state of the black hole 
are insensitive to the initial orbital eccentricity. 

C. Gravitational Waves 

1. General feature and convergence 

Figure [TT] plots gravitational waveforms as a function of a retarded time for runs APR1414H, APR1515H, 
APR1316H, and APR135165H. Here, the retarded time is defined by 

t,,t = t-D- 2Mo HD/Mo), (4.9) 
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TABLE III: Key numerical results for models APR145145, APR1515, APR1316, and APR135165. 
AE, A J, Afr>rAHi ^AH, Cp, Ce, aud /qnm denote energy and angular momentum carried by 
gravitational waves, rest mass of the material located outside apparent horizon, area of the apparent 
horizon in units of IGTrMgjjf, polar and equatorial circumferential radii of the apparent horizon, 
and quasinormal mode frequency, respectively. Mr>rAH ^^i^ /qnm are given in units of Mo and 
kHz, respectively, /qnm is derived from Eq. (|4.11|l using MeH.f and af2. For all the models, results 
for three grid resolutions are presented. All the quantities are evaluated when we stopped the 
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where D is the distance between the source and an observer. For comparison, inspiral waveforms calculated by the 
Taylor T4 formula are plotted together by the dashed curves. 

For the case that a hypermassive neutron star is the outcome after the merger (for model APR1414), gravita- 
tional waves are composed of the inspiral, merger, and quasiperiodic waveforms. Here, the merger waveform denotes 
short-term burst-type waves emitted after the inspiral phase and before the hypermassive neutron star relaxes to a 
quasisteady state. The quasiperiodic waveforms seen for model APR1414 are emitted by rapidly rotational motion 
of the ellipsoidal hypermassive neutron star and their amplitude is comparable to that of the late inspiral phase. 
Gravitational waves of nearly identical frequency (/ ~ 3.8 kHz for this model) are emitted for many cycles in the 
quasiperiodic phase, because the dissipation time scale of the rotational kinetic energy by gravitational radiation 
reaction is much longer than the rotational time scale. Consequently, the effective (time-integrated) amplitude of 
gravitational waves is much larger than that for the late inspiral phase of frequency / ~ 1 kHz (cf. Sec. IIV C 41) . This 
property is essentially the same as that found in our previous works [ll, [l^ . 

For the case that a black hole is formed in a dynamical time scale after the onset of merger (for models APR145145, 
APR1515, APR1316, and APR135165), gravitational waves are composed of the inspiral, merger, and ringdown 
waveforms. Because a rotating black hole in a quasisteady state is the final outcome after the merger, the ringdown 
waveform is characterized by the fundamental quasinormal mode of the formed black hole. This is a universal 
qualitative feature that holds irrespective of the total mass and mass ratio of the binary neutron stars, but the 
amplitude and characteristic frequency of gravitational waves depend strongly on these parameters (see Sec. IIV C 3p . 
All these properties were not clarified in the previous work [1]. The merger waveform denotes short-term burst-type 
waves emitted after the inspiral phase and before the ringdown phase. The frequency of gravitational waves in the 
merger phase is slightly lower than the frequency of the black-hole quasinormal mode, but the amplitude is much 
larger than that of the ringdown waveform. The frequency and amplitude of the merger waveform are seen in the 
Fourier space in a better manner (see discussion of Sees. IIV C 31 and IIV C i)) . 

2. Inspiral gravitational waves 

Figure [TT] shows that the inspiral waveforms are fitted in part well with those derived by the Taylor T4 formula. 
For comparing the frequency of gravitational waves as a function of time, the orbital angular velocity derived from 
gravitational waves is plotted in Fig. [12] for models APR1414, APR1515, APR1316, and APR135165. Here, the 
angular velocity, Q, is computed from ^'4 by 
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In Fig. [T2I the solid, dashed, short-dashed, and dotted curves are numerical results with three grid resolutions and 
results by the Taylor T4 formula, respectively. In comparison, the time axis for the results by the Taylor T4 formula 
is shifted according to the prescription described in Sec. IIIDI 

Figure [T^ shows that with the lowest-resolution run, the numerical results agree approximately with that of the 
Taylor T4 formula only for a short time (2 ms < irct ^ 4 ms). The time span, during which the numerical results are 
fitted well by the Taylor T4 formula, is longer for the better grid resolutions. Thus, disagreement of the numerical 
waveforms with those by the Taylor T4 formula is simply due to the fact that the grid resolution is insufficient. 
The more specific reason is that angular momentum is numerically dissipated more for the lower grid resolution, as 
mentioned in Sec. IIVBI As a result, merger time (defined as the time when the merger sets in, and approximately 
equal to the time spent in the inspiral phase) is spuriously shortened and the number of the wave cycles in the inspiral 
phase is spuriously reduced, resulting in a disagreement of the orbital phase with the Taylor-T4's result. These 
properties are observed in all the models (see also Fig. [0)1 . 

The numerical results for the angular velocity curve do not agree with the curve derived by the Taylor T4 formula 
even with the highest resolution run. To clarify the reason for this disagreement, we generate Fig. [T4j in which the 
results (a) for runs APR1515L and APR1515M and (b) for runs APR1316L and APR1316M are replotted by shifting 
the time axis to align the merger time for the results of three grid resolutions. Then, we also replot the result by 
the Taylor T4 formula so that it fits well with the overlapped numerical results. This figure shows that the results of 
three grid resolutions approximately agree each other and also with the result obtained by the Taylor T4 formula for 
4 ms < ij-et ^ 8 ms, whereas for t^ot ^ 8 ms three numerical results deviate from the result of the Taylor T4 formula. 
This suggests that (i) for the late inspiral phase (for 4 ms < t^-ct < 10 ms in this figure), the numerical results are 
approximately convergent and reliable and (ii) the disagreement with the results by the Taylor T4 formula comes 
from a reason different from the numerical error. The most plausible reason is that for such a late inspiral phase, a 
finite size effect of neutron stars, which is not taken into account in the Taylor T4 formula, plays an important role 
for accelerating the inward motion [63 | because tidal deformation of two neutron stars and its effect are not negligible 
for the last orbits. 

In binary BH system, the comparison of the PN and numerical-relativity waveforms have been extensively done, 
varying the PN order, changing the expression of the PN formula (e.g., Taylor Tl, T4, and Et), and including spins 
of black holes [H, [t^, [8l|, [12, HI] ■ The consensus is that the Taylor T4 formula is valid for mofl < 0.1. However, this 
is not the case for binary neutron stars because the finite size effect plays an important role for them, as mentioned 
above. We infer that the criterion will depend on EOS of neutron stars, because it determines compactness of neutron 
stars and thus degree of tidal deformation. For proving this point more strictly, however, it is necessary to compute 
waveforms with a higher precision in which a sufhcient convergence can be achieved. This is the issue left for the 
future. 

Before closing this subsection, we touch on the convergence of the inspiral waveforms. As mentioned above, the wave 
cycle is spuriously shorten for the lower grid- resolution simulations. Figure [15] (a) plots the merger time as a function 
of Ax^, where we define the merger time as the time when MqO reaches 0.12 (see Fig. [T2t . We find that the merger 
time converges approximately at second order. Comparing the merger time derived by the highest-resolution run 
with that derived by the extrapolated result, we find that the merger time is by '--^ 2 ms underestimated irrespective 
of models. This suggests that ^ 2 wave cycles are spuriously lost: Figure 13 (b) suggests that the wave cycles in 
an early phase of the numerical simulation seem to be lost. To derive more accurate gravitational waveforms for 
the inspiral phase, say, those with the error of half cycle, a higher grid resolution is required. For this purpose, the 
present simulation with uni-grid domain is computationally expensive (although it is in principle possible to perform 
the simulation), and hence, an adaptive mesh refinement (AMR) algorithm would be required [13, [H,!!^. Because 
the purpose of this paper is not to derive highly precise inspiral gravitational waveforms but to qualitatively explore 
gravitational waves in the merger and subsequent phases, we do not pay attention to the inspiral waveforms in more 
detail. 



3. Merger and ringdown gravitational waves 

In contrast to the results for the inspiral waveforms, the numerical results for the merger and ringdown waveforms 
have a good convergence, and the grid resolution chosen in this work appears to be acceptable. Figure fTSlfb) plots ^'4 
as a function of the retard time for runs APR1515L, APR1515M, and APR1515II. To focus on the ringdown waveforms 
in comparison among three runs, the time is shifted for runs APR1515L and APR1515M, to align the phase in the 
ringdown waveforms. This figure shows that the phase of the waveforms agrees well for 10 ms < t-^ct ^5 H-S ms. The 
amplitude of gravitational waves slightly disagree, but the error is at most ~ 20% for the lowest resolution run and for 
the highest resolution run, it appears to be less than 5% (note that this does not originate from the finite extraction 
radius as shown in Fig. [1]). This indicates that the grid resolution for model APR1515II and resulting spurious short 
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merger time do not seriously affect the merger and ringdown waveforms. This feature does not depend on the total 
mass or mass ratio of the model. Thus for a quantitative study of the merger and ringdown waveform, the highest 
grid resolution adopted in this paper is acceptable. 

As mentioned in Sec. IIV C 1[ the waveforms in the ringdown phase are primarily characterized by the fundamental 
quasinormal mode of the formed black holes for models APR145145, APR1515, APR1316, and APR135165. To clarify 
this fact, Fig. [16] plots \E'4 together with a fitting formula in the form 

Ae-*/'^sm{2nfQNMt + 6), (4.11) 

where A and 6 are constants, and the frequency and damping time scale are predicted by a linear perturbation analysis 

as [el [el 



/qnm « 10.7 (j^) [1 - 0.63(1 - flf)" '] kHz, 



(4.12) 



...^^^ms. (4.13) 

Figure [TBI fa) and (b) show that gravitational waves for model APR1515 are well fitted by the hypothetical waveforms 
given by Eq. (|4.1ip for ^^ct ^ 10.5 ms. [Here, we set Of = af2 for the fitting (cf. Sec. IIV B]) ]. This is reasonable because 
the final outcome for model APR1515 is a stationary rotating black hole with negligible disk mass, and hence, the 
black-hole perturbation theory (i.e., Eq. (|4.1ip ) should work well. This is also the case for the gravitational waveforms 
for model APR145145, for which the merger proceeds in essentially the same manner as that for model APR1515. 

We note that gravitational waves emitted for t^-ct < 10.5 ms are not well fitted by the hypothetical fitting formula. 
This implies that they are not the ringdown waves but merger waves. Namely, these gravitational waves are not 
emitted by the ringdown oscillation of the black hole but probably by a motion of the material moving around the 
black hole. Figure [TH] (a) and (b) show that the gravitational- wave amplitude in the merger phase is much larger than 
that in the ringdown phase. 

The fitting between the numerical and analytic waveforms for models APR1316 and APR135165 works fairly well, 
but is not as good as that for models APR1515 and APR145145 (see Fig. [T6| (c) and (d) for the results of run 
APR1316H); the numerical waveform may be fitted by the analytic one for tj-d ^ 10.4 ms, but the damping time 
appears to modulate. This disagreement is reasonable because for models APR1316 and APR135165, a fraction of 
material is located outside the black hole horizon even after formation of the black hole, and it subsequently falls into 
the black hole. Thus, the system is not completely in vacuum nor in a stationary state, and hence, the numerical 
waveforms may not be well fitted by the analytic results derived in an ideal assumption. 

Damping time scale td in Eq. (|4.13p allows us to evaluate black hole spin. By this method, we can check the accuracy 
of the spin parameters estimated by the methods described in Sec IIVBI Figure [iTj depicts the evolution of | ^'4 1 for 
runs APR1515 and APR1316, which shows a clear exponential decay in 10.5 ms ^ t < 11.5 ms for APR1515 and 
10 ms 10.5 ms for APR1316. We find that the damping time scale is estimated from this figure as 0.19 ± 0.01 

ms for run APR1515H and 0.20 ± 0.02 ms for run APR1316H. The larger uncertainty for run APR1316H reflects the 
fact that surrounding matter falls into black hole even after the black hole formation for a relatively long time scale 
and ringdown gravitational waves are not simply characterized by the fundamental quasinormal modes. Putting the 
values of the damping time scale into Eq. (|4.13p . the black hole spin is estimated as 0.77 ± 0.04 for run APR1515H 
and 0.84 ±0.06 for run APR1316H. As expected, these values agree with the results in Sec. IIV B] within the error bar. 

The merger waveforms also depend on the mass ratio. We plot the plus mode of gravitational waves (/i+) for 
runs APR1515H and APR135165H in Fig. [THl To focus on waveforms just before/after BH formation, waveforms 
arc shifted to align the formation time of apparent horizon (vertical dashed line). The peak at t ~ 10 ms for model 
APR1515 is higher than that for model APR135165. This reflects the fact that the less massive neutron star for model 
APR135165 is tidally deformed from the late inspiral phase, and at the merger, the material of it starts expanding. 
Consequently, its compactness decreases quickly, so does the degree of nonaxisymmetry of the system, reducing the 
amplitude of gravitational waves. 

The amplitude of the first peak after the black hole formation is, by contrast, higher for model APR135165 than 
for model APR1515. This seems to reflect the fact that more material is located outside the black-hole horizon at 
the black hole formation for model APR135165, because such material subsequently falls into the black hole to excite 
gravitational waves. 
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4- Fourier spectrum 



We define a Fourier power spectrum of gravitational waves by 



Hf) ^ / 



\h+ifW + \hA.fW 



2 



(4.14) 



where 



h+if) = j e 
hAf)= I e 



f'h+{t)dt, 



(4.15) 



(4.16) 



and /i+ and /ix denote the + and x modes of gravitational waves of I = \m\ = 2. Then, from h{f), we define a 
nondimensional spectrum (or effective amplitude) as 



Figure [H] shows the spectrum (hcs) of gravitational waves for runs APR1414H, APR1515H, APR1316H, and 
APR13516H. To plot Fig. [191 we assume D = 100 Mpc. Because the simulations are started with an orbit of a 
finite value of Hq {f = Hq/tt ^ 700 Hz), the spectrum amplitude is not realistic for a low- frequency side / < 800 
Hz. To compensate for this drawback, we plot the spectrum of gravitational waves derived by the Taylor T4 formula 
by the dotted curve, which approximately behaves as oc /~" where n depends weakly on / and is slightly larger 
than 1/6 around / < 1 kHz (note that for / ^ 1 kHz, n ^ 1/6 The figures show that the spectrum computed 
numerically smoothly connects to that of the Taylor T4 formula at ~ 0.8-1 kHz, indicating an acceptable accuracy 
of the numerical results. 

The first noteworthy feature found from Fig. [19] is that the spectrum amplitude does not steeply damp even at 
/ ^ 1.2-1.3 kHz j35] where the binary neutron stars reach the ISCO, although this is predicted by the results derived 
from the Taylor T4 formula. This indicates that gravitational waves are emitted by a high-velocity motion of a 
merging object even after the binary neutron stars reach the ISCO. Rather, the spectrum amplitude steeply damps 
at a fairly high frequency / = /cut ~ 2.5-3 kHz for the case that a black hole is formed. This indicates that by 
an inspiral-type motion, gravitational waves are emitted even up to such a high frequency: Even after the onset of 
merger, two high-density peaks remain in the merging object (see the fourth panel of Fig. 4), and the system emits 
gravitational waves for which the waveform is similar to the inspiral one. 

For the case that a hypermassive neutron star is formed (for model APR1414), multiple characteristic peaks for 
2-5 kHz are seen. These peaks are related to quasiperiodic gravitational waves emitted by the hypermassive neutron 
star. As reported in Ref. the peak of the highest amplitude appears at f ^ 3.8 kHz for model APR1414, which 
is associated with gravitational waves emitted by the quasiperiodic rotation of the hypermassive neutron star. The 
side-band peaks are generated by the coupling modulation between the modes of the quasiperiodic rotation and of 
a quasiradial oscillation of the hypermassive neutron star for which the oscillation frequency is ~ 1.2 kHz. We note 
that the simulation was artificially stopped at t ~ 20 ms for this model, but the hypermassive neutron star is likely 
to survive for a longer time > 1 s. Hence, the peak amplitude of the spectrum may be a factor of ~ 7('^ -^1000/20) 
larger in reality. In Fig. [20] (b), we plot h{f) together with the planned noise level of the advanced LIGO, to clarify 
detectability of quasiperiodic gravitational waves emitted by the hypermassive neutron star. Based on the hypothesis 
that the hypermassive neutron star survives for 1 s, the peak amplitudes around 3-4 kHz may be raised up by a 
factor of several. This suggests that such gravitational waves may be detected for an event of D < 20-30 Mpc, as 
mentioned in Refs. (ll. [23j. 

For models APR1515, APR145145, APR1316, and APR135165, the spectrum shape has qualitatively an universal 
feature (see Fig. [10] (a)): (i) for / < /cut ~ 2.5-3 kHz, the spectrum amplitude gradually decreases with / according 
to the relation cx /~" where n takes values between 1/6 and ^ 1/3; (ii) for / > /cut, the amplitude steeply decreases. 
This is likely to reflect the fact that two high-density peaks in the merged object disappear during the collapse to a 
black hole; (iii) for / = /peak ~ 5-6 kHz, which is slightly smaller than /qnm ~ 6.5-6.9 kHz, a broad peak appears. 
Because the frequency is always smaller than /qnm, this peak is not associated with the ringdown gravitational 
waveform but with the merger waveform; it is likely to be emitted by matter moving around the black hole; (iv) for 
/ > /fin ~ /qnm, the amplitude damps in an exponential manner. 

We note that the feature of the spectrum shape is qualitatively similar to that for the merger of binary black holes 
psj . except for the following differences: One of the most noticeable differences is found in the phase (ii). For the 
case of binary neutron stars, the spectrum amplitude steeply decreases for / > /cut- By contrast for binary black 



h^sif) = h{f)f. 



(4.17) 
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holes such a steep decrease is not found. The other remarkable difference is found in the peak amplitude associated 
with the phase (iii) and (iv): For the binary neutron stars, this is at most half of the amplitude at the ISCO (/ ~ 1 
kHz for the binary neutron stars), whereas for binary black holes, this peak amplitude is as large as or even larger 
than the amplitude at the ISCO (e.g., Ref. [is*!). 

The spectrum shape in the case of black hole formation is significantly different from that of hypcrmassive neutron 
star formation for / > 2 kHz. Figure [20l (b) compares the spectra for models APR1414 and APR145145. For model 
APR1414, there are multiple spiky peaks in the spectrum for 2 kHz ^ / < 6 kHz. By contrast, the spectrum shape 
is fairly smooth for model APR145145. This suggests that if gravitational waves of frequency between about 2 kHz 
and 6 kHz are detected, the outcome (hypermassive neutron star or black hole) can be distinguished. 

The most remarkable difference among the spectrum shapes of the four black-hole formation models shown in Fig. 
1201 (a) is found in the amplitude and width of peak associated with the merger waveform: The peak amplitude for 
models APR1515 and APR145145 is much less prominent than that for models APR135165 and APR1316 (see Fig. HOI 
(a) and compare the same-mass models), reflecting the fact that the amplitude of the merger gravitational waveform is 
smaller. This difference seems to reflect the difference in the disk-formation process, as mentioned in Sec. IIV C 31 For 
the equal-mass models, nearly all the material simultaneously collapse to a black hole, whereas for the unequal-mass 
models, disks surrounding the black hole are formed and subsequently the material of the disk fall into the black hole, 
enhancing the amplitude of gravitational waves. This difference in the accretion process is also reflected in the width 
of the peak: For the unequal-mass models, the width of the peak is broader. This indicates that a variety of the 
matter motion induces those gravitational waves. 

The merger process is also reflected in the value of /cut- For equal-mass runs APR145145H and APR1515H, 
/cut ~ 3.0 kHz, while for unequal-mass runs APR1316H and APR135165H, f^ut ~ 2.5 kHz (/cut does not depend 
strongly on the grid resolution). Thus, /cut is smaller for the unequal-mass models for a given value of the total mass. 
The reason for this is that tidal elongation and disruption of the less massive neutron star occur for the unequal-mass 
binary neutron stars in the late inspiral phase. The tidal elongation sets in at a relatively low frequency just after 
the binary neutron stars reach the ISCO. This is reflected in a steep decrease of the spectrum amplitude at a smaller 
frequency. 

The total mass of the system is reflected in the value of /gn. Figure HD] (a) shows that the frequency at which hctf 
reaches 2 x 10^^'^ (for D — 100 Mpc) is approximately inversely proportional to the total mass; e.g., /fin ^ 7.2 kHz 
for models APR1515 and APR135165 and /fi„ ~ 7.6 kHz for models APR145145 and APR1316. This is reasonable 
because the value of /fin is determined primarily by /qnm which is inversely proportional to the black-hole mass for 
a given spin, and hence, approximately to the total mass. 

Before closing this section, we comment on the convergence. Figure [20] (c) and (d) plots the spectrum for runs 
APR1515L, APR1515M, and APR1515H and that for runs APR135165L, APR135165M, and APR135165H, respec- 
tively. This shows that these do not overlap completely, but the amplitude for a given frequency agrees among three 
results within ~ 1 x 10"^'^ error for both models. Such error does not affect the findings and conclusions in this 
section. 



V. SUMMARY 

This paper reports new numerical results of general relativistic simulations for binary neutron stars, focusing in 
particular on the case that a black hole is formed in dynamical time scale after the onset of merger. Following Ref. [l|, 
the APR EOS and irrotational velocity field are employed for modeling the binary neutron stars in the inspiral phase. 
We prepare initial conditions with an orbital separation which is much larger than that in the previous work [Ij , and 
hence, unphysical effects associated with initial nonzero eccentricity and incorrect approaching velocity are excluded 
with a much better manner. We adopt the moving puncture approach, which enables to simulate black hole formation 
and subsequent longterm evolution of the black hole spacetime. Furthermore, the simulations are performed with a 
much better grid resolution than that in the previous work ! l] to obtain more reliable numerical results. 

As a result of the improvements summarized above, we have obtained the following new results for the outcome 
formed after the merger in the present work: (i) For the binary neutron stars modeled by the APR EOS, a black hole 
is formed in the dynamical time scale after the onset of merger, if the total mass of the system mp is larger than a 
threshold mass Mthr — 2.8-2.9Mq. This holds irrespective of the mass ratio as long as 0.8 < Qm £ 1- In the present 
work, the value of Afthr is determined with less uncertainty than that in Ref. [l| . (ii) For the case that the black hole 
is formed in the dynamical time scale after the onset of merger, the resulting black-hole spin Of is w 0.78 ±0.02. This 
value depends weakly on the total mass and mass ratio of the binary neutron stars, (iii) The mass of the formed black 
hole, MeHjf , is calculated from the approximate energy conservation relation (|4.3p and the circumferential equatorial 
radius of apparent horizon within 0.5% error, (iv) The mass and spin of the formed black hole are consistent with 
those determined by the frequency and damping time of the quasinormal mode of gravitational waves, even when a 
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disk of mass < O.OSM© is formed around the black hole, (v) A quasisteady disk of mass > O.OIMq is formed around 
the black hole for the mass ratio Qm k, 0.8. The disk mass depends not only on the mass ratio but also on the total 
mass of the binary neutron stars. We find that for the total mass (mo) closer to Mthr, the resulting disk mass is 
larger. 

Gravitational waves emitted during the merger and ringdown phases are also analyzed. For the case that a hyper- 
massive neutron star is formed, gravitational waves are composed of the inspiral, merger, and quasiperiodic waveforms. 
The quasiperiodic waves are emitted due to quasiperiodic rotation of the formed ellipsoidal hypermassive neutron star. 
For model APR1414, its characteristic frequency is ~ 3.8 kHz, which agrees with that found in the previous paper 
Q . As mentioned in the previous works ^ [2^ , the effective amplitude for the quasiperiodic waves could be larger 
than ~ 10"^"'^ at a distance of D < 100 Mpc. Although the frequency is rather high and slightly outside the best- 
sensitive band of ground-based gravitational- wave detectors, such gravitational waves are an interesting target for the 
next-generation interferometric detectors. 

For the case that a black hole is formed in the dynamical time scale after the onset of merger, gravitational waves 
are composed of the inspiral, merger, and ringdown waveforms. The feature is qualitatively universal, irrespective 
of the total mass and mass ratio of the binary neutron stars. The feature is clearly seen in the Fourier spectrum, 
and quantitatively summarized as follows: (i) For / < /cut ~ 2.5-3 kHz, the spectrum amplitude gradually decreases 
according to the relation cx where n is a slowly varying function of /; n = 1/6 for / <C 1 kHz, and n ~ 1/3 
for / — > /cut- /cut is much larger than the frequency at the ISCO. This is due to the fact that even after the onset 
of merger, the merging object has two high-density peaks, and emits gravitational waves for which the waveform is 
similar to the inspiral one. (ii) For / > /cut, the spectrum amplitude steeply decreases. This seems to reflect the fact 
that at such frequency, two density peaks disappear during the collapse to a black hole, (iii) For / = /peak ~ 5-6 kHz, 
which is slightly smaller than /qnm ~ 6.5-6.9 kHz, a broad peak appears. Because the frequency is always smaller 
than /qnm, this peak is not associated with the ringdown gravitational waveform but with the merger waveform, (iv) 
For / > /fin « /qnm, the amplitude damps in an exponential manner. 

Although the features (i)-(iv) are qualitatively universal, the values of /cut, /peak, /qnm, and /fin, and the height 
and width of the peak at f — /peak depend on the total mass and mass ratio of binary neutron stars, i.e., merger and 
black hole formation processes. This implies that if gravitational waves of high frequency / = 2-8 kHz are detected, 
we will be able to get information about merger and black hole formation precesses. 

As mentioned in Sec. I, the hybrid EOS adopted in this work seems to be appropriate for studying black hole 
formation in the merger of the binary neutron stars. However, for studying hypermassive neutron star formation 
or evolution of accretion disk around the formed black hole, the present choice of the EOS is not very appropriate, 
because for such cases, effects associated with the thermal energy (finite temperature) and neutrino cooling are likely 
to play an important role (e.g., Ref. [l^l). Magnetic fields will also play an important role for longterm evolution of 
the hypermassive neutron star and black hole accretion disks if they are amplified in a short time scale (say in 10 
dynamical time scales of the system). In the future paper, we would like to report our new studies in which these 
effects are taken into account. 

In the present work, we are not able to compute the inspiral waveforms with a high precision because the grid 
resolution is not high enough. Even with the highest-grid resolution (e.g., in run APR1515H), the phase error may 
amount to 2 wave cycles (see Fig. I13p . To suppress the phase error within, e.g., a half wave cycle in the inspiral phase, 
the grid resolution would have to be by 30-40% as high as that for model APR1515H. To perform such high-resolution 
simulation, the present choice of the grid structure requires a high computational cost, in particular for a systematic 
survey of gravitational waveforms for a variety of parameters. We have to adopt an AMR algorithm 0, [66| to 
save the computational cost. Currently, we are studying the inspiral waveforms in the simulation using a code with 
the AMR algorithm [s^l- We hope to present such results in the near future. 
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FIG. 1: The gravitational wave amplitude as a function of retarded time for run APR145145H. The solid, dashed, and dotted 
curves denote the amplitudes extracted at 513.4, 546.7, and 579.9 km, respectively. 




FIG. 2: The evolution of a coordinate separation defined by Eq. (|3.5|l for model APR1515. The solid, dashed, and short-dashed 
curves show the results for runs APR1515L, M, and H with approaching velocity, respectively. The dotted curve does for run 
APR1515L without approaching velocity. 
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FIG. 7: Snapshots of the density contour curves for p and velocity field («", w^) on the a; = plane for runs (a) APR1414H, (b) 
APR145145H, (c) APR1316H, and (d) APR135165H. The time shown in the upper-left side denotes the elapsed time from the 
beginning of the simulation. The short-dashed half circles on the equatorial plane denotes the location of the apparent horizon. 
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curves denote the results for low, medium, and high-resolution runs, respectively. 




FIG. 14: (a) The same as Fig. [IHb) but the numerical results for runs APR1515L and APR1515M are plotted by shifting the 
time coordinates to align the merger time with run APR1515H. (b) The same as (b) but for model APR1316. 




FIG. 15: (a) The numerical results for merger time as a function of As^ for all the models and (b) ^'4 as a function of retarded 
time for runs APR1515L, APR1515M, and APR1515H. To align the phase, the time coordinates for runs APR1515L and 
APR1515M are shifted. 
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FIG. 16: Ringdown waveforms associated primarily with the fundamental quasinormal mode (a), (b) for run APR1515H and 
(c), (d) for run APR1316H. The panels (a) and (c) plot the real part of *I'4 and (b) and (d) the imaginary part. For all 
the panels, the short-dashed curves denote the fitting curves calculated by Eq. (|4.1ip . The vertical dashed line in each panel 
denotes the formation time of the apparent horizon. 
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FIG. 17; Decay of 'I'4 for (a) run APR1515H and (b) run APR1316H. Solid, dashed, and short-dashed curves denote the real, 
imaginary part of ^'4, and exponential decay with td ~ 0.19 ms for run APR1515H and 0.20 ms for run APR1316H. 
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FIG. 18: The plus mode of gravitational waves {h+) for runs APR1515H (solid curve) and APR135165H (dashed curve). The 
vertical dashed line denotes the formation time of apparent horizon for run APR135165H. The waveform for run APR1515H is 
shifted to align the formation time of the apparent horizon. 
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FIG. 19: Spectrum of gravitational waves as a function of frequency for runs (a) APR1414H, (b) APR1515H, (c) APR1316H, 
and (d) APR135165H for a hypothetical distance of 100 Mpc. For comparison, the spectrum calculated by the Taylor T4 
(Newtonian quadrupole, e.g., oc f~^^^) formula is shown by the dashed curve (dotted line). 
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FIG. 20: (a) Spectrum of gravitational waves for runs APR145145H, APR1515H, APR1316H, and APR135165H, (b) spectrum 
for runs APR1414H and APR145145H, (c) the spectrum for runs APR1515L, APR1515M, and APR1515H, and (d) the spectrum 
for runs APR135165L, APR135165M, and APR135165H. Lines and cross symbols in the panel (a) denote /"^''^, f'^''^, and 
QNM frequencies for all the models. The dotted line in the panel (b) shows the planned noise level of the advanced LIGO. 



